{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d732c49a-eede-4ff4-8b04-08a4972df47e",
   "metadata": {},
   "source": [
    "# 作业3——快速傅里叶变换fft"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8fc686e7-945a-4a24-a329-36f94471dd99",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAASq0lEQVR4nO3dbYxcZ3mH8etm7YDFi9zKdDGbgKlqRQpEqmEVYkWVRpTWiRs1FqJS+ECAVjKhQSIqBAhIIPVDebFEAYXasgBBVBpAYEyETF2qZARINpA4L8YJWwwU4heRApqEhRXEy90Pc+zdTGY9s97ZnLPPXj9p5JnnPDNzz73P/nNy5sxOZCaSpJXvGXUXIEkaDQNdkgphoEtSIQx0SSqEgS5JhVhT1xNv2LAhN23aVNfTn/Ob3/yGZz/72XWX0Qj2Yo69mGMv5jShF/fee+8vMvP5/bbVFuibNm3innvuqevpz2m327RarbrLaAR7McdezLEXc5rQi4j46ULbPOQiSYUw0CWpEAa6JBXCQJekQhjoklSIgWe5RMSzgG8Cz6zmfykz398zJ4CPAduB3wJvzMwjoy9Xpdp/30l2HZziVGeGF65fxy3bLmXHlom6y6qFvdCFGua0xd8Br8rM6YhYC3w7Ir6emYfnzbkG2FxdXgnsrv6VBtp/30lu3XeUmSdmATjZmeHWfUcBVl2Q2QstxcBDLtk1Xd1cW116/+budcDt1dzDwPqI2DjaUlWqXQenzgXYWTNPzLLr4FRNFdXHXmgphvpgUUSMAfcCfwZ8IjO/0zNlAnhk3u0T1djpnsfZCewEGB8fp91uX1jVIzQ9Pd2IOpqgrl6c7MwsOF7Xz8ZezPF3ZE7TezFUoGfmLPDnEbEe+EpEvCwzvz9vSvS7W5/H2QvsBZicnMy6P3EFzfjkV1PU1YuJw3f1DbKJ9etq+9nYizn+jsxpei8WdZZLZnaANnB1z6YTwCXzbl8MnFpKYVo9btl2KevWjj1pbN3aMW7ZdmlNFdXHXmgpBgZ6RDy/2jMnItYBrwZ+0DPtTuCG6LoSeCwzTyMNYceWCT7wmsu5aKy7HCfWr+MDr7l8Vb4JaC+0FMMcctkIfLY6jv4M4IuZ+bWIuBEgM/cAB+iesnic7mmLb1qmelWoHVsmuOO7PwPgC2/eWnM19bIXulADAz0zHwS29BnfM+96AjeNtjRJ0mL4SVFJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhBgZ6RFwSEXdHxMMRcSwi3tZnTisiHouI+6vL+5anXEnSQtYMMecM8PbMPBIRzwXujYhvZOZDPfO+lZnXjr5ESdIwBu6hZ+bpzDxSXf818DAwsdyFSZIWZ5g99HMiYhOwBfhOn81bI+IB4BTwjsw81uf+O4GdAOPj47Tb7cXWO3LT09ONqKMJ6u5FpzMD0Iifh72YU3cvmqTpvRg60CPiOcCXgZsz8/GezUeAF2fmdERsB/YDm3sfIzP3AnsBJicns9VqXWDZo9Nut2lCHU1Qdy92Tx0CoNXaWlsNZ9mLOXX3okma3ouhznKJiLV0w/xzmbmvd3tmPp6Z09X1A8DaiNgw0kolSec1zFkuAXwKeDgzP7LAnBdU84iIK6rH/eUoC5Uknd8wh1yuAl4PHI2I+6ux9wAvAsjMPcBrgbdExBlgBrg+M3P05UqSFjIw0DPz20AMmHMbcNuoipIkLZ6fFJWkQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBViYKBHxCURcXdEPBwRxyLibX3mRER8PCKOR8SDEfHy5SlXy2H/fSe56oN38cb//A1XffAu9t93su6S1ACui5VnzRBzzgBvz8wjEfFc4N6I+EZmPjRvzjXA5urySmB39a8abv99J7l131FmnpgF4GRnhlv3HQVgx5aJOktTjVwXK9PAPfTMPJ2ZR6rrvwYeBnp/otcBt2fXYWB9RGwcebUauV0Hp8790p4188Qsuw5O1VSRmsB1sTINs4d+TkRsArYA3+nZNAE8Mu/2iWrsdM/9dwI7AcbHx2m324urdhlMT083oo66nOzMLDj+dPelU9XShJ9H3eui7l40aV00Sd3rYpChAz0ingN8Gbg5Mx/v3dznLvmUgcy9wF6AycnJbLVaw1e6TNrtNk2ooy4Th+/q+8s7sX7d096X3VOHAGi1tj6tz9tP3eui7l40aV00Sd3rYpChznKJiLV0w/xzmbmvz5QTwCXzbl8MnFp6eVput2y7lHVrx540tm7tGLdsu7SmitQErouVaZizXAL4FPBwZn5kgWl3AjdUZ7tcCTyWmacXmKsG2bFlgg+85nIuGusuhYn16/jAay73ja9VznWxMg1zyOUq4PXA0Yi4vxp7D/AigMzcAxwAtgPHgd8Cbxp5pVo2O7ZMcMd3f0an0+Hgu15VdzlqCNfFyjMw0DPz2/Q/Rj5/TgI3jaooSdLi+UlRSSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIdYMmhARnwauBR7NzJf12d4Cvgr8pBral5n/PMIapaHsv+8kuw5OcaozwwvXr+OWbZeyY8vEqqtBq9fAQAc+A9wG3H6eOd/KzGtHUpF0Afbfd5Jb9x1l5olZAE52Zrh131GApy1Qm1CDVreBh1wy85vAr56GWqQLtuvg1LkgPWvmiVl2HZxaVTVodRtmD30YWyPiAeAU8I7MPNZvUkTsBHYCjI+P0263R/T0F256eroRddSt05lhdna21l50OjMAF1TDyeq+/cYv5PEuZF2Msoal9GKUmrAumqTpeTGKQD8CvDgzpyNiO7Af2NxvYmbuBfYCTE5OZqvVGsHTL0273aYJddRt99QhOp1Orb3YPXUIgFZr66LvO3H4rr6BOrF+3QW9pgtZF6OsYSm9GKUmrIsmaXpeLPksl8x8PDOnq+sHgLURsWHJlUmLcMu2S1m3duxJY+vWjnHLtktXVQ1a3Za8hx4RLwB+npkZEVfQ/Y/EL5dcmbQIZ990fOeXHuT3s39gooYzTJpQg1a3YU5bvANoARsi4gTwfmAtQGbuAV4LvCUizgAzwPWZmctWsbSAHVsmuOO7PwPgC2+u51BFE2rQ6jUw0DPzdQO230b3tEZJUo38pKgkFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQAwM9Ij4dEY9GxPcX2B4R8fGIOB4RD0bEy0dfpiRpkGH20D8DXH2e7dcAm6vLTmD30suSJC3WwEDPzG8CvzrPlOuA27PrMLA+IjaOqkBJ0nDWjOAxJoBH5t0+UY2d7p0YETvp7sUzPj5Ou90ewdMvzfT0dCPqqFunM8Ps7Gytveh0ZgCWVMMoHgOWti6a9DqWqgnrokmanhejCPToM5b9JmbmXmAvwOTkZLZarRE8/dK0222aUEfddk8dotPp1NqL3VOHAGi1ttb6GLC0ddGk17FUTVgXTdL0vBjFWS4ngEvm3b4YODWCx5UkLcIoAv1O4IbqbJcrgccy8ymHWyRJy2vgIZeIuANoARsi4gTwfmAtQGbuAQ4A24HjwG+BNy1XsZKkhQ0M9Mx83YDtCdw0sookSRfET4pKUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKMVSgR8TVETEVEccj4t19trci4rGIuL+6vG/0pUqSzmfNoAkRMQZ8Avgr4ATwvYi4MzMf6pn6rcy8dhlqlCQNYZg99CuA45n548z8PfB54LrlLUuStFgD99CBCeCRebdPAK/sM29rRDwAnALekZnHeidExE5gJ8D4+DjtdnvRBY/a9PR0I+qoW6czw+zsbK296HRmAJZUwygeA5a2Lpr0OpaqCeuiSZqeF8MEevQZy57bR4AXZ+Z0RGwH9gObn3KnzL3AXoDJyclstVqLKnY5tNttmlBH3XZPHaLT6dTai91ThwBotbbW+hiwtHXRpNexVE1YF03S9LwY5pDLCeCSebcvprsXfk5mPp6Z09X1A8DaiNgwsiolSQMNE+jfAzZHxEsi4iLgeuDO+RMi4gUREdX1K6rH/eWoi5UkLWzgIZfMPBMRbwUOAmPApzPzWETcWG3fA7wWeEtEnAFmgOszs/ewjCRpGQ1zDP3sYZQDPWN75l2/DbhttKVJkhbDT4pKUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKEZk5eFLE1cDHgDHgk5n5wZ7tUW3fDvwWeGNmHjnfYz5z4+bc+IaPXmDZkrQ6nf7szfzu9A+j37aBe+gRMQZ8ArgGuAx4XURc1jPtGmBzddkJ7F5SxZKkRVszxJwrgOOZ+WOAiPg8cB3w0Lw51wG3Z3d3/3BErI+IjZl5eqEHvfjXj/Iv3/q3JZQuSavPDefZNswx9AngkXm3T1Rji51DROyMiHsi4p4hnleStAjD7KH3O1bTe+B9mDlk5l5gL3SPob/rL/5xiKeXJJ3z45sX3DTMHvoJ4JJ5ty8GTl3AHEnSMhom0L8HbI6Il0TERcD1wJ09c+4EboiuK4HHznf8XJK0DDJz4IXu6Yj/A/wIeG81diNwY3U96J4J8yPgKDA56DFf8YpXZBPcfffddZfQGPZijr2YYy/mNKEXwD25QK4OcwydzDwAHOgZ2zPvegI3LfG/LZKkJfCTopJUCANdkgphoEtSIQx0SSrEUH+ca1meOOL/gJ/W8uRPtgH4Rd1FNIS9mGMv5tiLOU3oxYsz8/n9NtQW6E0REfdk5mTddTSBvZhjL+bYizlN74WHXCSpEAa6JBXCQK/+WJgAezGfvZhjL+Y0uher/hi6JJXCPXRJKoSBLkmFWHWBHhHviIiMiA3zxm6NiOMRMRUR2+aNvyIijlbbPl59GfaKFxG7IuIHEfFgRHwlItbP27aqetErIq6uXvvxiHh33fUst4i4JCLujoiHI+JYRLytGv/jiPhGRPyw+veP5t2n7xopQUSMRcR9EfG16vbK6sNCf4axxAvdL+E4SPcDTRuqscuAB4BnAi+h+yeAx6pt3wW20v3zwF8Hrqn7NYyoD38NrKmufwj40GrtRU9fxqrX/KfARVUvLqu7rmV+zRuBl1fXn0v3z2RfBnwYeHc1/u5h1kgJF+CfgP8AvlbdXlF9WG176P8KvJMnfz3edcDnM/N3mfkT4DhwRURsBJ6XmYey+xO8HdjxdBe8HDLzvzLzTHXzMN1vmIJV2Ise574QPTN/D5z9QvRiZebpzDxSXf818DDd7wO+DvhsNe2zzP28+66Rp7XoZRIRFwN/A3xy3vCK6sOqCfSI+FvgZGY+0LNpoS+4nqiu946X5u/p7nGDvRjqy85LFRGbgC3Ad4DxrL51rPr3T6ppJffoo3R3+P4wb2xF9WGoL7hYKSLiv4EX9Nn0XuA9dA81POVufcbyPOMrwvl6kZlfrea8FzgDfO7s3frMX/G9WITV8jqfIiKeA3wZuDkzHz/PWyRF9igirgUezcx7I6I1zF36jNXeh6ICPTNf3W88Ii6ne5zrgWqhXgwciYgrWPgLrk8wdyhi/viKsFAvzoqINwDXAn9ZHUaBQnuxCKvyy84jYi3dMP9cZu6rhn8eERsz83R1yO3RarzUHl0F/G1EbAeeBTwvIv6dldaHug/i13EB/pe5N0VfypPf3Pgxc28Efg+4krk3ArfXXfuIXv/VwEPA83vGV10vel7/muo1v4S5N0VfWnddy/yag+57Ih/tGd/Fk98M/PCgNVLKBWgx96boiupDUXvoFyIzj0XEF+kG3BngpsycrTa/BfgMsI5uiH2974OsPLfRXYjfqP6P5XBm3rhKe3FOZp6JiLfSPRNqDPh0Zh6ruazldhXweuBoRNxfjb0H+CDwxYj4B+BnwN/BwN+XEq2oPvjRf0kqxKo5y0WSSmegS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEL8Pyon0S/Ub0HlAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 720x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#例\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "\n",
    "def getSin(amp,freq,phase,sampleList):\n",
    "    return amp*np.sin(2*math.pi*freq*sampleList+phase)\n",
    "\n",
    "# 1. 获得混合波形\n",
    "srate=100000\n",
    "t=np.linspace(0,1,srate)\n",
    "s1=getSin(amp=1.5, freq=50, phase=0,sampleList=t)    \n",
    "s2=getSin(amp=3, freq=100, phase=0,sampleList=t)  \n",
    "s3=getSin(amp=2, freq=200, phase=0,sampleList=t) \n",
    "m=s1+s2+s3\n",
    "\n",
    "# 2. 获得傅里叶系数\n",
    "fCoefs = np.fft.fft(m,srate)\n",
    "\n",
    "# 3. 获得振幅列表：每一个绕线的重心到原点的距离\n",
    "amp_list=2*np.abs(fCoefs/srate)\n",
    "\n",
    "# 把频率轴从0~1000 转变成 0~499 然后 -500~-1\n",
    "freqs = np.fft.fftfreq(len(amp_list), 1/srate)\n",
    "\n",
    "# 然后把 频率轴 和 数据 都变成 0hz 在中间，向左是负频率，向右是正频率的形式\n",
    "amp_shifted=np.fft.fftshift(amp_list)\n",
    "freq_shift=np.fft.fftshift(freqs)\n",
    "\n",
    "# 把振幅列表画出来\n",
    "plt.grid()\n",
    "plt.stem(freq_shift,amp_shifted) #注意：由于频率有负值，所以振幅被平分了，需要乘以2才能还原原始振幅\n",
    "plt.xlim([-500,500])\n",
    "plt.figure(figsize=(10,4))   \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "1de4be9b-25bb-4445-bedb-da1570dcafa4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAAD4CAYAAACg9uHUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAesklEQVR4nO3df5DkdX3n8ec7y5JbkbvVMBAYSMBzmQqJFzehIBSV1JDkBLZMWL3oQVIRL9at5mJVsJQEohVNUndyoczlUhjIqhR4MQiVwErixpUT5tASRGAXFsQNCxrc2T22UAcZmdNled8f853YO3Rvz/b32/P9fnuej6qu7v7+6O+rZz7f97y/3f3ticxEkiRJ9fihugNIkiStZDZjkiRJNbIZkyRJqpHNmCRJUo1sxiRJkmp0VN0Bulm7dm2++tWvrjtGX9/97nc55phj6o7Rlzmr15asbcn5wAMPPJOZY3XnqIL1q1ptyQntyWrO6pWtYY1sxk444QTuv//+umP0NTU1xeTkZN0x+jJn9dqStS05I+Kf685QFetXtdqSE9qT1ZzVK1vDfJtSkiSpRjZjkiRJNbIZkyRJqpHNmCRJUo1sxiRJkmrUtxmLiFMi4q6IeCwiHo2I3y2mvzIi7oiIx4vrV/RY/4KI2BURuyPiiqqfwJHasn2ac6+6k9Ou+DTnXnUnW7ZP1x1JaqRR2FdGrX5JGk1LeWXsBeDdmfkTwM8BvxMRZwBXAJ/LzHXA54r7h4iIVcCHgQuBM4BLinVrsWX7NFfeupPpmTkSmJ6Z48pbd7byj4w0TCO0r4xM/ZI0uvo2Y5m5LzMfLG4/BzwGjAMXATcWi90IbOyy+lnA7sx8MjO/D3yyWK8WV2/bxdyBg4dMmztwkKu37aopkdRMo7KvjFL9kjS6juhLXyPiVGA98CXghMzcB/MFLyKO77LKOPCNjvt7gLN7PPYmYBPA2NgYU1NTRxJtSaZn5npOH2R7s7OzQ8lZNXNWry1ZB81Z9b7SBG2vX1Ub9TFch7ZkNWfzLLkZi4iXA38HXJaZ34mIJa3WZVp2WzAzNwObASYmJnIY37o7fu+dXf/IjK9dM9C3/Lbl24HNWb22ZB00Z9X7St1GoX5VbdTHcB3aktWczbOksykjYjXzhewTmXlrMfnpiDixmH8isL/LqnuAUzrunwzsHTxuOZefP8Ga1asOmbZm9SouP3+ipkRSM43SvjIq9UvS6FrK2ZQBfAx4LDP/rGPW7cClxe1LgU91Wf3LwLqIOC0ijgYuLtarxcb143zwja/h6FXzT3t87Ro++MbXsHH9eF2RpEYalX1llOqXpNG1lLcpzwV+E9gZETuKaX8AXAXcEhFvA54C3gQQEScBH83MDZn5QkS8E9gGrAKuz8xHK34OR2Tj+nFuuu8pAG5++zl1RpEabUT2lZGqX5JGU99mLDO/QPfPTgD8Upfl9wIbOu5vBbYOGlCSBmX9ktQGfgO/JElSjWzGJEmSamQzJkmSVCObMUmSpBrZjEmSJNXIZkySJKlGNmOSJEk1shmTJEmqkc2YJElSjWzGJEmSamQzJkmSVCObMUmSpBrZjEmSJNXIZkySJKlGNmOSJEk1OqrfAhFxPfB6YH9m/lQx7WZgolhkLTCTma/tsu7XgeeAg8ALmXlmJaklSZJGRN9mDLgBuAb4+MKEzPyPC7cj4kPAs4dZ/7zMfGbQgJJUhgeUkpqubzOWmXdHxKnd5kVEAG8GfrHiXJJUlRvwgFJSg5X9zNjPA09n5uM95ifw2Yh4ICI2ldyWJB2xzLwb+Fa3eR0HlDctayhJ6rCUtykP5xIOX8TOzcy9EXE8cEdEfLUojC9RNGubAMbGxpiamioZrbeZmTmA0tuYnZ0das6qmLN6bclaNmdV+0qDLfWAMoG/yszN3RZazvpVlZUyhpdTW7Kas4Eys+8FOBV4ZNG0o4CngZOX+BgfAN6zlGVPP/30HKY3X/fFfPN1Xyz9OHfddVf5MMvAnNVrS9ayOavaV/oB7s8l1IZBL91qWDH9WuDdh1nvpOL6eOAh4Bf6bWvY9asqK2UML6e2ZDVn9crWsDJvU/4y8NXM3NNtZkQcExHHLtwGXgc8UmJ7klSZiDgKeCNwc69lMnNvcb0fuA04a3nSSVpJ+jZjEXETcA8wERF7IuJtxayLWfQWZUScFBFbi7snAF+IiIeA+4BPZ+ZnqosuSaV4QCmpEZZyNuUlPaa/tcu0vcCG4vaTwE+XzCdJpRQHlJPAcRGxB3h/Zn6MHgeUwEczcwPzB5S3zX/Gn6OAv/GAUtIwlP0AvyQ1mgeUkprOf4ckSZJUI5sxSZKkGtmMSZIk1chmTJIkqUY2Y5IkSTWyGZMkSaqRzZgkSVKNbMYkSZJqZDMmSZJUI5sxSZKkGtmMSZIk1chmTJIkqUY2Y5IkSTWyGZMkSaqRzZgkSVKN+jZjEXF9ROyPiEc6pn0gIqYjYkdx2dBj3QsiYldE7I6IK6oMLkmSNAqW8srYDcAFXab/j8x8bXHZunhmRKwCPgxcCJwBXBIRZ5QJK0lHygNKSU3XtxnLzLuBbw3w2GcBuzPzycz8PvBJ4KIBHkeSyrgBDyglNdhRJdZ9Z0S8BbgfeHdmfnvR/HHgGx339wBn93qwiNgEbAIYGxtjamqqRLTDm5mZAyi9jdnZ2aHmrIo5q9eWrGVzVrWv1Ckz746IUwdY9V8OKAEiYuGA8isVxpOkgZuxa4E/AbK4/hDwW4uWiS7rZa8HzMzNwGaAiYmJnJycHDBaf9fuugeAyclzSj3O1NQUw8xZFXNWry1Zy+asal9pqMoOKJfzYLIqK+WAYjm1Jas5m2egZiwzn164HREfAf6hy2J7gFM67p8M7B1ke5JUsUoPKJfzYLIqK+WAYjm1Jas5m2egr7aIiBM77r4BeKTLYl8G1kXEaRFxNHAxcPsg25OkKmXm05l5MDNfBD7C/FuSi3lAKWlZLOWrLW4C7gEmImJPRLwN+NOI2BkRDwPnAe8qlj0pIrYCZOYLwDuBbcBjwC2Z+eiQnockLZkHlJKapO/blJl5SZfJH+ux7F5gQ8f9rcBLzlKSpOVSHFBOAsdFxB7g/cBkRLyW+bcdvw68vVj2JOCjmbkhM1+IiIUDylXA9R5QShqGMmdTSlLjeUApqen8d0iSJEk1shmTJEmqkc2YJElSjWzGJEmSamQzJkmSVCObMUmSpBrZjEmSJNXIZkySJKlGNmOSJEk1shmTJEmqkc2YJElSjWzGJEmSamQzJkmSVCObMUmSpBrZjEmSJNWobzMWEddHxP6IeKRj2tUR8dWIeDgibouItT3W/XpE7IyIHRFxf4W5JWlJrGGSmm4pr4zdAFywaNodwE9l5r8D/gm48jDrn5eZr83MMweLKEml3IA1TFKD9W3GMvNu4FuLpn02M18o7t4LnDyEbJJUmjVMUtNFZvZfKOJU4B8y86e6zPt74ObM/Osu874GfBtI4K8yc/NhtrEJ2AQwNjb2s7fccstSn8MR++CX5gC48uw1pR5ndnaWl7/85VVEGipzVq8tWcvmrGpf6ee88857YJivPA27hi1n/arKShnDy6ktWc1ZvdI1LDP7XoBTgUe6TH8vcBtFU9dl/knF9fHAQ8AvLGV7p59+eg7Tm6/7Yr75ui+Wfpy77rqrfJhlYM7qtSVr2ZxV7Sv9APfnEmrDoJflrGHDrl9VWSljeDm1Jas5q1e2hg18NmVEXAq8HviNIki3Rm9vcb2/KHhnDbo9SaqSNUxSUwzUjEXEBcDvA7+amc/3WOaYiDh24TbwOuCRbstK0nKyhklqkqV8tcVNwD3ARETsiYi3AdcAxwJ3FKd8X1cse1JEbC1WPQH4QkQ8BNwHfDozPzOUZyFJPVjDJDXdUf0WyMxLukz+WI9l9wIbittPAj9dKp0klWQNk9R0fgO/JElSjWzGJEmSamQzJkmSVCObMUmSpBrZjEmSJNXIZkySJKlGNmOSJEk1shmTJEmqkc2YJElSjWzGJEmSamQzJkmSVCObMUmSpBrZjEmSJNXIZkySJKlGNmOSJEk1shmTJEmqUd9mLCKuj4j9EfFIx7RXRsQdEfF4cf2KHuteEBG7ImJ3RFxRZfA6bdk+zblX3clbP/Ndzr3qTrZsn647knQIx+gPWMMO5dhQ063EMbqUV8ZuAC5YNO0K4HOZuQ74XHH/EBGxCvgwcCFwBnBJRJxRKm0DbNk+zZW37mR6Zg6A6Zk5rrx154oYLGoHx+hL3IA1DHBsqPlW6hjt24xl5t3AtxZNvgi4sbh9I7Cxy6pnAbsz88nM/D7wyWK9Vrt62y7mDhw8ZNrcgYNcvW1XTYmkQzlGD2UN+wHHhppupY7RowZc74TM3AeQmfsi4vguy4wD3+i4vwc4u9cDRsQmYBPA2NgYU1NTA0brb6bouAfZxkK33m36MDOXMTs729hsndqSE5qdtcoxWmZfabhKa9hy1q8yrF/D1ZasTc7ZxjFahUGbsaWILtOy18KZuRnYDDAxMZGTk5NDigXX7roHgMnJc4543fF77+w6WMbXrmGYmcuYmppqbLZObckJzc5a5Rgts6+MgCXXsOWsX2VYv4arLVmbnLONY7QKg55N+XREnAhQXO/vsswe4JSO+ycDewfcXmNcfv4Ea1avOmTamtWruPz8iZoSSYdyjC7Jiqxhjg013Uodo4M2Y7cDlxa3LwU+1WWZLwPrIuK0iDgauLhYr9U2rh/ng298DUevmv/Rja9dwwff+Bo2rh+vOZk0zzG6JCuyhjk21HQrdYz2fZsyIm4CJoHjImIP8H7gKuCWiHgb8BTwpmLZk4CPZuaGzHwhIt4JbANWAddn5qPDeRrLa+P6cW667ylmZmbY9vu/WHcc6SUcoz9gDTuUY0NNtxLHaN9mLDMv6THrl7osuxfY0HF/K7B14HSSVJI1TFLT+Q38kiRJNbIZkyRJqpHNmCRJUo1sxiRJkmpkMyZJklQjmzFJkqQa2YxJkiTVyGZMkiSpRjZjkiRJNbIZkyRJqpHNmCRJUo1sxiRJkmpkMyZJklQjmzFJkqQa2YxJkiTVaOBmLCImImJHx+U7EXHZomUmI+LZjmX+sHRiSZKkEXLUoCtm5i7gtQARsQqYBm7rsujnM/P1g25HkoYhIiaAmzsmvQr4w8z8845lJoFPAV8rJt2amX+8TBElrRADN2OL/BLwRGb+c0WPJ0lD5QGlpKao6jNjFwM39Zh3TkQ8FBH/GBE/WdH2JKlKHlBKqk3pV8Yi4mjgV4Eru8x+EPjxzJyNiA3AFmBdj8fZBGwCGBsbY2pqqmy0nmZm5gBKbWNmZo6DBw8ONWdVZmdnzVmxNmStYoxWsa+0RN8DSmAv8J7MfHTxAstZv6pg/RqOtmRtQ842jdEqVPE25YXAg5n59OIZmfmdjttbI+IvI+K4zHymy7Kbgc0AExMTOTk5WUG07q7ddQ8Ak5PnlHqMmZkZhpmzKlNTU+asWBuyVjFGq9hXmq6KA8rlrF9VsH4NR1uytiFnm8ZoFap4m/ISehxRRsSPRkQUt88qtvfNCrYpSVU57AFlZs4Wt7cCqyPiuOUOKGm0lXplLCJeBvx74O0d094BkJnXAb8G/HZEvADMARdnZpbZpiRV7LAHlMDTmZkeUEoallLNWGY+D/zIomnXddy+BrimzDYkaVg8oJTUBFV9tYUktY4HlJKawH+HJEmSVCObMUmSpBrZjEmSJNXIZkySJKlGNmOSJEk1shmTJEmqkc2YJElSjWzGJEmSamQzJkmSVCObMUmSpBrZjEmSJNXIZkySJKlGNmOSJEk1shmTJEmqkc2YJElSjUo1YxHx9YjYGRE7IuL+LvMjIv4iInZHxMMR8TNltidJVbKGSWqCoyp4jPMy85ke8y4E1hWXs4Fri2tJagprmKRaDfttyouAj+e8e4G1EXHikLcpSVWxhkkaurKvjCXw2YhI4K8yc/Oi+ePANzru7ymm7Vv8QBGxCdgEMDY2xtTUVMlovc3MzAGU2sbMzBwHDx4cas6qzM7OmrNibchaxRitYl9puEpq2HLWrypYv4ajLVnbkLNNY7QKZZuxczNzb0QcD9wREV/NzLs75keXdbLbAxVFcDPAxMRETk5OlozW27W77gFgcvKcUo8xMzPDMHNWZWpqypwVa0PWKsZoFftKw1VSw5azflXB+jUcbcnahpxtGqNVKPU2ZWbuLa73A7cBZy1aZA9wSsf9k4G9ZbYpSVWxhklqgoGbsYg4JiKOXbgNvA54ZNFitwNvKc5I+jng2cx8yVuUkrTcrGGSmqLM25QnALdFxMLj/E1mfiYi3gGQmdcBW4ENwG7geeA/lYsrSZWxhklqhIGbscx8EvjpLtOv67idwO8Mug1JGhZrmKSm8Bv4JUmSamQzJkmSVCObMUmSpBrZjEmSJNXIZkySJKlGNmOSJEk1shmTJEmqkc2YJElSjWzGJEmSamQzJkmSVCObMUmSpBqV+Ufh0kC2bJ/m6m272Dszx0lr13D5+RNsXD9+xMtIUh361Sfrl46UzZiW1Zbt01x5607mDhwEYHpmjitv3QnA2iUsY0GTVKd+Ncz6pUH4NqWW1dXbdv1LkVowd+AgV2/bdUTLSFId+tUn65cGYTOmZbV3Zq7v9KUsI0l16FefrF8ahM2YltVJa9f0nb6UZSSpDv3qk/VLgxi4GYuIUyLiroh4LCIejYjf7bLMZEQ8GxE7issflourtrv8/AnWrF51yLQ1q1dx+fkTR7SMVJY1TIPoV5+sXxpEmQ/wvwC8OzMfjIhjgQci4o7M/Mqi5T6fma8vsR2NkIUPsP7e3z7M9w++yHjHmUZTU4/3XUaqkDVMR6xfDbN+aRADN2OZuQ/YV9x+LiIeA8aBxYVMOsTG9ePcdN9TANz89nMGXkYqwxqmQfWrT9YvHalKvtoiIk4F1gNf6jL7nIh4CNgLvCczH+3xGJuATQBjY2NMTU1VEa2rmeKDlGW2MTMzx8GDB4easyqzs7ONy9ntd7A4ZxW/p2Fp4s90sSrGaJN/B1UqW8OWs35VwfpVXr8a1uR9p6k/005tGqNVKN2MRcTLgb8DLsvM7yya/SDw45k5GxEbgC3Aum6Pk5mbgc0AExMTOTk5WTZaT9fuugeAycnBj1iu3XUPMzMzDDNnVaamphqXs9vvYHHOKn5Pw9LEn+liVYzRJv8OqlJFDVvO+lUF61d5/WpYk/edpv5MO7VpjFah1NmUEbGa+SL2icy8dfH8zPxOZs4Wt7cCqyPiuDLblKSqWMMkNUGZsykD+BjwWGb+WY9lfrRYjog4q9jeNwfdpiRVxRomqSnKvE15LvCbwM6I2FFM+wPgxwAy8zrg14DfjogXgDng4szMEtuUpKpYwyQ1QpmzKb8ARJ9lrgGuGXQbkjQs1jBJTeE38EuSJNWokq+20Ojasn2aq7ftYu/MHCe1/MsLR+m5SOpvlPb5UXoueimbMfW0Zfs0V966k7kDBwGYnpnjylt3ArSuCIzSc5HU3yjt86P0XNSdb1Oqp6u37fqXnX/B3IGDXL1tV02JBjdKz0VSf6O0z4/Sc1F3NmPqaW/xDdJLnd5ko/RcJPU3Svv8KD0XdWczpp5OWrvmiKY32Sg9F0n9jdI+P0rPRd3ZjKmny8+fYM3qVYdMW7N6FZefP1FTosGN0nOR1N8o7fOj9FzUnR/gV08LHwz9vb99mO8ffJHxFp/BM0rPRVJ/o7TPj9JzUXcrqhlbODV4unif/dQrPu2g7mPj+nFuuu8pAG5+e/P+4e2RGKXnMmzuKxoFo7TPj9Jz0UutmGZs8anBCzxFWDqU+4okLa8V85mxbqcGL/AUYekH3FckaXmtmGas3ynAniIszXNfkaTltWKasbUvW33Y+f9mzeHnL9iyfZpzr7qTL33tWzzx7Its2T5dRTypcs889z2eePZFTrvi05x71Z1LHqtV7SuSpKVZEc3Y+7bs5NvPHzjsMt/5fwf6/rFa+CzNwoeaX3gRLrt5B+v/+LM2ZWqMLdunee0ffZYnnvkuL7wIyQ8+79VvnFa1r0iSlm7km7H3bdnJX9/7VN/lXkx41807DvtH5o/+/tGun6X59vMHeNfNO3jflp2lskplvW/LTt518w5m5l7aUM0dOMgHbn+057q/8ZF7KttXJElL18qzKedfoXqYuQMvVvq4yfwrXZfdvGOgdf/63qde8sfsmKNX8V/f8BrPPlNltmyf5gO3P9q14epnZu4Ap17x6dIZjmRfecXLVvP+X/lJ9wFJ6iEyc/CVIy4A/iewCvhoZl61aH4U8zcAzwNvzcwH+z3uD5+4Lk+89M8HziWpXfbdeBnf2/d4LPd2h1HDrF/SylO2hg38NmVErAI+DFwInAFcEhFnLFrsQmBdcdkEXDvo9iSpStYwSU1R5m3Ks4DdmfkkQER8ErgI+ErHMhcBH8/5l9/ujYi1EXFiZu473AOf/Nx+/tvn/7JENElt8pZ6NjuUGmb9klaesjWszAf4x4FvdNzfU0w70mUAiIhNEXF/RNxfItOSvWz1sr8jIrXOy1YHa//VyO4rldWw5a5fkkZLmVfGulXoxR9AW8oy8xMzNwObYf4zF7//8/+lRLTeFj5Q/yvrx4d2IoA0Cs79t6/kE/95/n/gDX1fefKy4Tzu4VVWw5arfklqqJI1rEwztgc4peP+ycDeAZZZFt3O6Nq4fpyNRVN2JGenLTR0X3nsK9zyeA50Vps0TD8U8Otn/xgve/7/HvEYrXJfabhW1TBJo6tMM/ZlYF1EnAZMAxcDv75omduBdxafxTgbeLbf58WORFWnzC/8oTlSa599nD/49cmu80bsj5Yaqt8+MDX1zZ5jdBBHsq+0YB+ovYZJEgCZOfCF+dO9/wl4AnhvMe0dwDuK28H82UpPADuBM5fyuKeffnq2wV133VV3hCUxZ/XakrUtOYH7s0QtGvQyjBpm/apWW3JmtierOatXtoaV+tLXzNwKbF007bqO2wn8TpltSNKwWMMkNcHI/zskSZKkJrMZkyRJqpHNmCRJUo1sxiRJkmpU6h+FD0tEPAfsqjvHEhwHPFN3iCUwZ/XakrUtOScy89i6Q1TB+lW5tuSE9mQ1Z/VK1bBSZ1MO0a7MPLPuEP1ExP3mrE5bckJ7srYpZ90ZKmT9qlBbckJ7spqzemVrmG9TSpIk1chmTJIkqUZNbcY21x1gicxZrbbkhPZkNefya8tzMWf12pLVnNUrlbWRH+CXJElaKZr6ypgkSdKKYDMmSZJUo1qbsYh4U0Q8GhEvRsSZHdNPjYi5iNhRXK7rmPezEbEzInZHxF9ERNSVs5h3ZZFlV0ScX2fOLrk/EBHTHT/HDf1y1yUiLiiy7I6IK+rO0ykivl78LncsnL4cEa+MiDsi4vHi+hU1Zbs+IvZHxCMd03pmq+v33iNna8ZnN22pX4fLWsxrZA1r0/iwfg2czfq1IDNruwA/AUwAU8CZHdNPBR7psc59wDlAAP8IXFhjzjOAh4AfBk4DngBW1ZWzS+4PAO/pMr1n7prGwaoiw6uAo4tsZ9Q5Nhfl+zpw3KJpfwpcUdy+AvjvNWX7BeBnOveXXtnq/L33yNmK8XmY59SK+tUna2NrWFvGh/WrVDbrV3Gp9ZWxzHwsM5f8TdURcSLwrzPznpx/1h8HNg4r34LD5LwI+GRmfi8zvwbsBs6qK+cR6Jq7xjxnAbsz88nM/D7wySJjk10E3FjcvpGafr+ZeTfwrUWTe2Wr7ffeI2cvTRufXbWlfsHI1bCmjQ/r14CsXz/Q5M+MnRYR2yPi/0TEzxfTxoE9HcvsKabVZRz4Rsf9hTxNyvnOiHi4eJl14eXeXrnr0rQ8iyXw2Yh4ICI2FdNOyMx9AMX18bWle6le2Zr4c27D+BxEG+oXNL+GtWF8NC3PYtav4alsfA793yFFxP8GfrTLrPdm5qd6rLYP+LHM/GZE/CywJSJ+kvmXyxer5Ls5BszZK8/Qcr4kwGFyA9cCf1Js+0+ADwG/tZz5lqhpeRY7NzP3RsTxwB0R8dW6Aw2oaT/nxo/PttQvaGcNs34tC+vXcFQ6PofejGXmLw+wzveA7xW3H4iIJ4DTme8wT+5Y9GRgb105izyndMkztJyLLTV3RHwE+Ifibq/cdWlankNk5t7ien9E3Mb8S85PR8SJmbmveEtnf60hD9UrW6N+zpn59MLtpo7PttSvYlutq2HWr+Gzfg1H1fWrkW9TRsRYRKwqbr8KWAc8Wbxk+VxE/FxxZs9bgF5HfMvhduDiiPjhiDityHlfU3IWA3nBG4CFM0G65l7ufB2+DKyLiNMi4mjg4iJj7SLimIg4duE28Drmf463A5cWi11KveNwsV7ZGvV7b9H4PCItql/Q4BrWovFh/arWyqxfy3EmQq9L8QT2MH8U+TSwrZj+H4BHmT8j4UHgVzrWObN40k8A11D8F4E6chbz3ltk2UXH2UZ15OyS+38BO4GHiwFyYr/cNY6FDcA/FZneW3eejlyvKsbhQ8WYfG8x/UeAzwGPF9evrCnfTcy/LXagGKNvO1y2un7vPXK2Znz2eE6tqF+Hy3q4n3XdNaxN48P6NXA+61dx8d8hSZIk1aiRb1NKkiStFDZjkiRJNbIZkyRJqpHNmCRJUo1sxiRJkmpkMyZJklQjmzFJkqQa/X+7+gYNl9JBugAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "\n",
    "# 定义一个函数，用于生成正弦波\n",
    "def getSin(amp, freq, phase, sampleList):\n",
    "    return amp * np.sin(-2 * math.pi * freq * sampleList + phase)\n",
    "\n",
    "# 定义一个函数，用于生成余弦波\n",
    "def getCos(amp, freq, phase, sampleList):\n",
    "    return amp * np.cos(-2 * math.pi * freq * sampleList + phase) \n",
    "\n",
    "# 定义去噪函数，将振幅小于阈值的噪声去除\n",
    "def denoise(arr, thresh):\n",
    "    mask = arr > thresh\n",
    "    return mask * arr\n",
    "\n",
    "# 1. 获得混合波形\n",
    "srate=3000\n",
    "t=np.linspace(0,1,srate)\n",
    "s1 = getSin(amp=1.5, freq=30, phase=0, sampleList=t)\n",
    "s2 = getCos(amp=3, freq=5, phase=0, sampleList=t)\n",
    "s3 = getSin(amp=10, freq=100, phase=0, sampleList=t)\n",
    "s4 = getCos(amp=20, freq=120, phase=0, sampleList=t)\n",
    "# 将四个正弦和余弦波相加，得到混合信号\n",
    "m = s1 + s2 + s3 + s4\n",
    "# 2. 获得傅里叶系数\n",
    "fCoefs = np.fft.fft(m,srate)\n",
    "\n",
    "# 3. 获得振幅列表：每一个绕线的重心到原点的距离\n",
    "amp_list=2 * np.abs(fCoefs / srate)\n",
    "\n",
    "# 把频率轴从0~1000 转变成 0~499 然后 -500~-1\n",
    "freqs = np.fft.fftfreq(len(amp_list), 1/srate)\n",
    "# 然后把 频率轴 和 数据 都变成 0hz 在中间，向左是负频率，向右是正频率的形式\n",
    "amp_shifted=np.fft.fftshift(amp_list)\n",
    "freq_shift=np.fft.fftshift(freqs)\n",
    "# 创建两个子图，分别用于显示原始振幅频谱和处理后的振幅频谱\n",
    "fig, ax = plt.subplots(1, 2, figsize=(10, 4))\n",
    "# 绘制原始振幅频谱\n",
    "ax[0].stem(freq_shift, amp_shifted)\n",
    "ax[0].set_xlim([-150, 150])\n",
    "#ax[0].set_ylim([-1, 21])\n",
    "ax[0].grid()\n",
    "# 对振幅频谱进行去噪处理\n",
    "# 1. 去掉所有低振幅噪音\n",
    "# mask = np.ones_like(amp_shifted) #创建一个全1数组作为掩膜\n",
    "# mask[np.abs(freq_shift) < 1] = 0 #将频率小于1的掩膜位置设为0\n",
    "# amp_shifted *= mask   #去除频率小于1的所有噪音\n",
    "amp_shifted = denoise(amp_shifted, 1)\n",
    "# 2. 去掉所有的高于110hz的正负频率噪音\n",
    "amp_shifted[(freq_shift > 110) | (freq_shift < -110)] = 0\n",
    "# 绘制处理后的振幅频谱\n",
    "ax[1].stem(freq_shift, amp_shifted)\n",
    "# 设置第二个子图的坐标轴范围\n",
    "ax[1].set_xlim([-150, 150])\n",
    "ax[1].set_ylim([-1, 21])\n",
    "# 给第二个子图添加网格线\n",
    "ax[1].grid()\n",
    "# 显示图形\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
